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Abstract. Several studies have shown that sub-mesoscales 
(SM ~ 1 km horizontal scale) play an important role in mixed 
layer dynamics. In particular, high resolution simulations 
have shown that in the case of strong down-front wind, the 
re-stratification induced by the SM is of the same order of the 
de-stratification induced by small scale turbulence, as well as 
of that induced by the Ekman velocity. These studies have 
further concluded that it has become necessary to include 
SM in ocean global circulation models (OGCMs), especially 
those used in climate studies. 

The goal of our work is to derive and assess an analytic 
parameterization of the vertical tracer flux under baroclinic 
instabilities and wind of arbitrary directions and strength. To 
achieve this goal, we have divided the problem into two parts: 
first, in this work we derive and assess a parameterization of 
the SM vertical flux of an arbitrary tracer for ocean codes that 
resolve mesoscales, M, but not sub-mesoscales, SM. In Part 
2, presented elsewhere, we have used the results of this work 
to derive a parameterization of SM fluxes for ocean codes 
that do not resolve either M or SM. 

To carry out the first part of our work, we solve the SM dy- 
namic equations including the non-linear terms for which we 
employ a closure developed and assessed in previous work. 
We present a detailed analysis for down-front and up-front 
winds with the following results: 

(a) down-front wind (blowing in the direction of the sur- 
face geostrophic velocity) is the most favorable condition for 
generating vigorous SM eddies; the de-stratifying effect of 
the mean flow and re-stratifying effect of SM almost cancel 
each other out. 
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(b) in the up-front wind case (blowing in the direction op- 
posite to the surface geostrophic velocity), strong winds pre- 
vents the SM generation while weak winds hinder the pro- 
cess but the eddies amplify the re-stratifying effect of the 
mean velocity, 

(c) wind orthogonal to the geostrophic velocity. In this 
case, which was not considered in numerical simulations, we 
show that when the wind direction coincides with that of the 
horizontal buoyancy gradient, SM eddies are generated and 
their re-stratifying effect partly cancels the de-stratifying ef- 
fect of the mean velocity. The case when wind direction is 
opposite to that of the horizontal buoyancy gradient, is anal- 
ogous to the case of up-front winds. 

In conclusion, the new multifaceted implications on the 
mixed layer stratification caused by the interplay of both 
strength and directions of the wind in relation to the buoy- 
ancy gradient disclosed by high resolution simulations have 
been reproduced by the present model. 

The present results can be used in OGCMs that resolve M 
but not SM. 


1 Introduction 

Recently, there has been a considerable interest in sub- 
mesoscales (SM) which are oceanic structures with sizes 
0(1 km) and a life time of the order of days. If one considers 
that the highest resolution 0(1/10°) in stand-alone OGCMs 
(ocean global circulation models) can represent structures of 
about 10 km which is 10 times larger than SM sizes and that 
OGCMs employed in thousand years runs for climate studies 
have in general a 1° resolution (corresponding to structures 
100 times larger than SM sizes), it seems clear that a good 
deal of important physical processes have thus far gone un- 
represented in many OGCMs. 

The parameterization of SM cannot be constructed by 
analogy with mesoscales, M. Indeed, while mesoscales are 
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characterized by a small Rossby number Ro=t; //<5C 1 (where 
£ and / are the relative and planetary vorticities respec- 
tively), mixed layer SM are characterized by 1 (Thomas 
et al., 2008). Thus, the dynamics of M and SM are quite dif- 
ferent and so are the final results for the fluxes of interest. For 
this reason, the parameterizations for M that were suggested 
for the deep ocean cannot be extrapolated to describe SM in 
the mixed layer and a new parameterization is required. 

Most of the present knowledge about mixed layer SM 
comes from high resolution numerical simulations (Levy et 
al., 2001; Thomas and Lee, 2005; Mahadevan, 2006; Ma- 
hadevan and Tandon, 2006; Klein et al., 2008; Thomas et al., 
2008; Capet et al., 2008, C8; Levy et al., 2010; Mahadevan 
et al., 2010, MTF). These studies have revealed many inter- 
esting features of SM, e.g., their contribution to the vertical 
mixing of buoyancy and tracers in the upper ocean. Among 
the most salient effects of SM on global ocean properties is 
the well documented tendency to re-stratify the mixed layer 
(Spall, 1995; Nurser and Zhang, 2000; C8; MTF). The effect 
of SM on deep convection has been recently demonstrated 
by Levy et al. (2009, Fig. 9) who point out a better agree- 
ment with the new mixed layer data by Boyer Montegut et 
al. (2004) south of the WBC (Western Boundary Current). 
Even in non-convective regimes, one can expect a signifi- 
cant cancellation between SM and small scale fluxes lead- 
ing to the mixed layer re-stratification (e.g., C8, Fig. 12; 
Klein et al., 2008); Hosegood et al. (2008) estimated that 
SM contribute up to 40% of the re-stratification process. 
In addition, as noted by Lapeyre et al. (2006) and Klein et 
al. (2008), the surface layers re-stratification is compensated 
by de-stratification of the ocean interior pointing to an in- 
teresting dynamical connection between surface and interior 
processes. Another important effect of SM concerns the lo- 
cation of the WBC that is shifted south by 4° and whose 
off shore extension penetrates further to the east, in better 
agreement with observations (Levy et al., 2009). An earlier 
study by Treguier et al. (2005), who found a significant in- 
crease (from ~30 Sv to ~70 Sv) in the barotropic transport 
in the Gulf Stream when moving from 1° to 1/6° resolution, 
was recently confirmed and the inclusion of SM further in- 
creased the transport by ~50 Sv. Finally, the structure of the 
MOC (meridional overturning circulation) was also signifi- 
cantly affected by SM not so much in its intensity as to its 
location (Levy et al., 2009, Fig. 12). 

This brief summary of some of the results of very high 
resolution (1/54°, ~2km) regional studies highlights the im- 
portance of SM and thus the question arises as to how much 
of SM physics is actually accounted for by global OGCMs. 
Even today’s highest resolution global ocean models ~1/10° 
(Maltrud and McClean, 2005; Sasaki et al., 2008) are not 
able to capture the SM field and much less are in the position 
to do so the OGCMs coupled to an atmospheric model used 
in climate studies where the resolution is 0(1°) but gener- 
ally lower. The significant global processes revealed in going 
from 1° resolution (~100km) and 1/54° resolution (~2km) 


are presently absent in such global models especially in cli- 
mate studies. Therefore, a reliable parameterization of SM in 
terms of the resolved fields has become necessary to ensure 
the physical completeness of mixed layer mixing processes. 

As MTF stressed, “parameterization of the circulation in- 
duced by SM eddies in the presence of wind forcing is 
required in climate models in order to simulate the re- 
stratification correctly”. To accomplish such a goal, a pa- 
rameterization: (1) must be valid for an arbitrary tracer since 
a model for buoyancy only is insufficient because it cannot 
describe an important ingredient such as CO 2 , (2) must in- 
clude a wind of arbitrary direction and intensity since MTF 
have shown that its effect on SM fluxes is large and because 
forcing in future climates is likely to be quite different from 
today’s, (3) must reproduce simulation data and finally, (4) to 
be usable in climate codes, it must be expressed only in terms 
of resolved fields, that is, the parameterization must be aver- 
aged over mesoscale fields. No parameterization presently 
available satisfies these criteria. 

We begin by considering the model independent dynamic 
equation for an arbitrary mean tracer 1 : 

At + Vtf • F h + d z F v = d z (K v d z r) + G (la) 

where the SM horizontal and vertical tracer fluxes are defined 
as follows: 

Fu = u'r' , Fy = w'r' (lb) 

In Eqs. (la,b), an overbar indicates averages over subme- 
soscales. To derive a parameterization for OGCMs that do 
not resolve either SM or M, one must further average the 
present results over the mesoscale fields, a problem we have 
studied elsewhere (Canuto and Dubovikov, 2010). 

The method we employ to derive the SM parameteriza- 
tion is analytical and thus it can be followed and checked in 
detail. The final result for the vertical SM flux for an arbi- 
trary tracer under arbitrary buoyancy and wind conditions, is 
also expressed analytically. The model includes non-linear 
interactions. To obtain the desired results, one carries out 
three steps: first, one solves in Fourier space the SM dy- 
namic equations describing the i/', w' , x' fields; second, one 
constructs the second-order correlation function w'x' in k- 
space and third, one integrates the results over all wave vec- 
tors to obtain the fluxes in physical space in terms of the 
resolved fields, to be used in Eq. (la). The procedure was 
first worked out for the linear case by Killworth (1997) and 
for the non-linear case by Canuto and Dubovikov (2005, 
2006, CD5, 6). Though the dynamic equations describing 
the SM velocity and temperature fields are formally the same 
as those describing mixed layer mesoscales that were dis- 
cussed in Canuto et al. (2010), in the present case they must 

1 A=3r+«-V ff+uJd z where (m, w) represent the mean flow, 
small scale vertical mixing is represented by the first term in the 
right hand side of Eq. (la) where K v is the vertical diffusivity, G 
represents external sources and primes denote submesoscale fields. 
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be solved in the regime appropriate to SM, namely for a 
Rossby number Ro=f /f= 0(1) rather than Ro-^l, as in the 
case of mesoscales. In addition, SM are trapped in the ML 
while mesoscales extend throughout the entire water column 
and form coherent structures (Provenzale, 1999). 

The key difficulty in solving the SM dynamic equations is 
the presence of non-linear terms whose closure is expressed 
in Eq. (3a) below. Since the latter is a key ingredient of the 
present model and since the original derivation (Canuto and 
Dubovikov, 1997) is somewhat involved, in Appendices A 
and B we have attempted to find a way to present a more 
physical approach to Eq. (3a) with the goal of highlighting 
the physical rather than the technical features of Eq. (3a). 

In addition to the derivation of the closure relations 
(Eq. 3a), there is the issue of the assessment of Eq. (3a) 
when applied to flows different than the present one so as 
to justify their use in the present context. Such an assess- 
ment was made using data from freely decaying flows, 2-D 
flows, rotating flows, unstably stratified flows, shear driven 
flows, DNS data, etc. and the results were in good agree- 
ment with the data (see Canuto et ah, 1999 and references 
therein). Even so, we consider such an assessment necessary 
but not sufficient for the credibility of the parameterization of 
SM fluxes derived below. The additional requirement con- 
sists in assessing the model predictions against results from 
SM resolving simulations. A first simulation corresponds to 
a system forced only by baroclinic instabilities and no wind 
(Fox-Kemper et ah, 2008, FFH) while a second one corre- 
sponds to a flow under realistic wind and buoyancy forcing 
(C8, 0.75 km resolution; MTF, 1 km resolution). They will 
be discussed in Sects. 5-7. 

The following two conditions must be further satisfied by 
an SM parameterization: (a) it must reproduce existing data 
and (b) it must predict new features to be assessed when such 
data become available. In this context, it must be mentioned 
that our work was posted as an OS Discussions (Canuto and 
Dubovikov, 2009, CD9) before Dr. A. Mahadevan kindly 
sent us the MTF manuscript. Our model would have been 
falsified had its predictions turned out to be inconsistent with 
MTF data. However, the model predictions in CD9 not only 
did not contradict the simulation data, but called attention to 
the same qualitative SM effects as MTF did in their paper. 

To make the SM parameterization usable in OGCMs, we 
looked for analytical solutions of the SM dynamic equations 
and to achieve that goal, we introduced the assumption that 
the fluxes are mostly contributed by their spectra in the vicin- 
ity of their maxima. Though this introduces errors of several 
tens of a percent, the advantages of obtaining analytic expres- 
sions for the vertical tracer flux in terms of resolved fields in 
the presence of both frontogenesis and Ekman pumping, was 
worth exploring. Another approximation which has helped 
us obtain analytical results follows from the assumption that 
the SM kinetic energy A/sm exceeds K—ii 2 / 2 where u is the 
baroclinic component of the mean velocity (we call atten- 
tion to the fact that K is considerably smaller than the mean 


kinetic energy), that is, the condition of applicability of the 
present treatment is predicated on the assumption: 

K <5C /fsM do) 

which we shall check several times in the following. Fur- 
thermore, following Killworth (2005), we adopt the approx- 
imation that due to the mixed layers strong mixing, one can 
neglect x z in the SM equations. Anticipating our main result, 
the vertical flux that enters Eq. (la) will be shown to have the 
following form: 

d z F y = h+ ■ V#t (Id) 

where uf plays the role of a bolus velocity. Since r- is small, 
one may make the analogy with the mesoscale bolus velocity 
more complete by adding to Eq. (Id) the term w$x z , where 
Wg is found from the continuity condition 3 z w s+^H *«§" =0 
(Killworth, 2005). 

The organization of the paper is as follows. In sec. 2 we 
discuss the dynamic equations for the SM fields in the ML 
and apply the turbulence closure model to the non-linear 
terms; in Sect. 3 we present the form of 9- Fy and Fy which 
we derive in Appendix C; in Sect. 4 we derive the explicit 
form of the SM kinetic energy K$m in terms of the resolved 
fields that, together with the results of the previous section, 
completes the problem of expressing 3 z Fy in terms of re- 
solved fields in the presence of both frontogenesis and Ek- 
man pumping. In Sect. 5 we study the case of a strong wind 
when the Ekman velocity exceeds the geostrophic one. We 
shall show that when a strong wind blows in the direction 
of the geostrophic velocity or of V#£>, it tends to de-stratify 
the mixed layer but at the same time it generates SM that 
tend to re-stratify the mixed layer. On the other hand, when 
the wind blows in directions opposite to geostrophic velocity 
or to V# fo, it re-stratifies the mixed layer, an effect that is 
strengthened by the re-stratifying effect of SM. In Sect. 6 we 
compare the model results with the data from the SM resolv- 
ing simulations of Capet et al. (2008). In Sect. 7 we compare 
the model results for the no-wind case. In Sect. 8, we present 
some conclusions. 

2 Sub-mesoscales dynamic equations near the surface 

Consider an arbitrary tracer field r and separate it into mean 
and fluctuating partsr=r+r / . The dynamical equation for 
the SM tracer field x' is obtained by subtracting the equation 
for the mean tracer r from that of the total field r. Since 
this procedure is well known and entails only algebraic steps 
with no physical assumptions, we cite only the final result 
(the notation is explained in footnote 1): 

At' = -V • Vr - Gh - Qv + dz (Xvdzx') (2a) 

Qh = « ' • r ' - «' ■ VtfT', Qy = w'x' z - w' x' z 

where the functions Q’s represent the non-linear terms. As 
expected, the average of Eqs. (2a) yields identically zero. It 
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must be noted that in Eqs. (2a) no closure has been used 
for the non-linear terms. Without the non-linear terms, 
Eqs. (2a) formally coincide with those describing mixed 
layer mesoscales tracer fields studied by Killworth (2005). 
The difference in representing M and SM lies in the scales 
over which the averages (represented by an overbar in 
Eqs. 2a), is taken: in the case of mesoscales, averages are 
over scales exceeding mesoscales while in the case of subme- 
soscales, averages are meant to be over scales smaller than 
mesoscales but larger than submesoscales. Furthermore, in 
describing mesoscales, one has Ro-^l and Ri'^> 1 , whereas in 
the case of sub-mesoscales, both Ro, Ri~0( 1). Following 
Killworth (2005), we neglect the terms containing r, and r', 
in which case the first of Eqs. (2a) simplifies to: 

d,r' + u ■ V H r' = - u ■ S7 H r - Q ^ (2b) 

Without the non-linear term, Eq. (2b) coincides with Eq. (2) 
of Killworth (2005) for the mesoscale buoyancy field. Within 
the same approximation, the equation for the horizontal SM 
velocity is given by: 

d t u' +u-TJ H u' +u' -V H u+ f e z xu' — -p _1 V# p'-Qn (2c) 

<2h = u ■ V#m' — u’ ■ V#h' (2d) 

where e z is the unit vector along z axis. Next, we Fourier 
transform Eqs. (2b,c) in horizontal planes and time. Follow- 
ing Killworth (1997, 2005), we keep the same notation u', r' 
for the submesoscale fields in the k— co space and assume that 
when Eqs. (2b,c) are Fourier transformed, the mean fields ti 
and Vtf r are constant in time and horizontal coordinates. We 
thus obtain: 

i{k ■ u — co)r — — u! ■ V#r — < 2 h 

i(k ■ u — co)ii = — u! ■ Vhu — fe z xu ’ — Q ^ — ikp~ l p’ 
d z w' = — Vfj ■ u = — ik ■ u (2e) 


where we have added the continuity equation that provides 
the z-derivative of w'. We recall that r', u’ and the non- 
linear terms are functions of the horizontal wave vector and 
frequency (k, co) and z, while u is a function of z only and 
V h r is z independent. 

Equations (2e) form a closed system whose solution pro- 
vides the necessary ingredients to construct the vertical flux 
(Eq. lb) provided one has a closure for the non-linear terms, 
a problem discussed in Appendices A and B with the result 
that, in the vicinity of k — ko where the SM energy spectrum 
E(k) has its maximum, the non-linear terms Q\\ have the 
following forms: 

Q T H (k, co) = xt \ k, co), Qn(k, co) = /;/(£, co), (3a) 


X — koUsM, 


K SM = 



where the scale £=kg 1 may be interpreted as the SM hori- 
zontal length scale. As it was shown in detail in CD5, ko is 
obtained from the solution of the eigenvalue problem which 
is derived from the eddy dynamical equations (Eq. 2e). In 
the limit of a strong non-linearity represented by: 

K sm >K, K = Ur (3b) 

where K is the kinetic energy of the baroclinic component 
of the mean velocity (defined in Eq. 4b), the solution of the 
eigenvalue problem yields the following result: 

k~ l =l^rs = n-\N/\f\)h (3c) 

where r$ is the Rossby deformation radius of the mixed layer 
(ML) of depth h and where N is the buoyancy frequency in 
the ML. Relation (Eq. 3c) is in qualitative agreement with 
other evaluations of the submesoscale length scale discussed 
in the literature (Boccaletti et ah, 2007; Thomas et ah, 2008; 
Fox-Kemper and Ferrari, 2008). In Fig. 9 of Fox-Kemper 
et ah (2008), the authors, using simulation data, plot SM 
length scales defined using different variables and in unit of 
the Stone length scale which is of the order of the deforma- 
tion radius. 

Equation (2e) together with Eq. (3a), represent a stochastic 
Langevin equation which has played a major role in turbu- 
lence modeling studies (Kraichnan, 1971; Leith, 1971; Her- 
ring and Kraichnan, 1971; Chasnov, 1991). The advantage 
of the Langevin equation is that it is linear in the fluctuat- 
ing fields and thus allows one to compute second-order mo- 
ments while the original Eqs. (2b,c) are non-linear and do 
not allow an analytical computation of such correlation func- 
tions. The key problem is to find a model for the non-linear 
terms Q’s that leads to a Langevin equation whose corre- 
lation functions are sufficiently close to those of the orig- 
inal Eqs. (2b,c). This is the closure problem for the non- 
linear terms. In CD5, we used the closure (Eq. 3a) derived 
by Canuto and Dubovikov (1997) and solved the eigenvalue 
problem to which the mesoscale dynamic equations were 
shown to reduce. Closure (Eq. 3a) has a simple interpreta- 
tion within the mixing length approach. In fact, the first two 
relations are quite standard with / _1 being the characteristic 
time scale while the third relation containing the character- 
istic length scale (Eq. 3c) and velocity, is the only possible 
combination that leads to a time scale. 

The advantage of the Langevin equation is that it allows 
us to express all SM fields in terms of the SM horizontal ve- 
locity u'(k,<o) which, in turn, allows us to express the spec- 
trum of any second-order moment in terms of the SM energy 
spectrum and of the resolved fields. Such a program, which 
we discuss in detail in Appendix C, was previously used by 
the authors to parameterize mesoscales in the ocean interior 
(CD5, 6) and in the mixed layer (Canuto et ah, 2010). The 
main difference in the case of SM is in the treatment of the 
eddy velocity field. In particular, the rotational component 
of the SM horizontal velocity called i/r in Eq. (C9), does 
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not coincide with the geostrophic velocity. In fact, Eq. (C12) 
shows that, using the third of Eq. (C6), in the limit (Eq. lc), 
hr can be represented as follows: 


Hr = - 


ikp 1 p' 

1 +\Ro 2 ' 


Ro = 


w 


(3d) 


where Ro is the Rossby number and p’ is the SM pressure 
field. Thus, when treating mesoscales that are character- 
ized by small Ro, the first relation in Eq. (3d) reduces to 
the geostrophic relation hr— >M g =— ikp~ l p’ . On other hand, 
since SM are characterized by Ro~l, we must use the com- 
plete form of Eq. (C12). This is the reason why we have 
not called hr the geostrophic component and the divergent 
component mr> (Eq. Cll) the a-geostrophic one. 


section. Before doing so and for future reference, we next 
derive the explicit form of the vertical flux itself. Integrating 
Eq. (4a) over z with the boundary condition Fy{ 0)=0, we ac- 
count for only the z-dependency of u within the mixed layer. 
Thus, we obtain: 


Fy — — kh ■ V#r 


(6a) 


where the submesoscale diffusivity is given by: 
kh = z(l + K 2 ) 

Z 

udz 
o 



u — y ~ — x u 

’ I/I 


(6b) 


3 Sub-mesoscale vertical tracer flux 

Following the program we have outlined at the end of the 
previous section, in Appendix C we show how to express 
w'(k, &>) and r'{ k, o>) in terms of the SM horizontal velocity 
u\k , ft;) and of the resolved fields; that, in turn, allows us to 
express the spectrum of Fy = w'r' in terms of the SM energy 
spectrum. Finally, integrating the spectra over all wavenum- 
bers, we obtain the following SM vertical tracer flux 

3 rFv=Wg -V//T, =— ^l+y 2 ^ 

where: 

rsl/l 1 ~ _ if 

Y = = — , u — u — h I u(z)dz, (4b) 

UsM Ro J 

—h 

e z is the unit vertical vector and Ro is the Rossby number 
defined in Eq. (3d). It is worth stressing that in the second 
relation in Eq. (4a), the second term in the square bracket 
is a vector: in fact, although c z xm is a pseudo-vector (cross 
product of the vectors e z and a), f is a pseudo-scalar which 
is the scalar product of the vector e z and the pseudo-vector 
2 £2 and thus the product is a vector. The variable li may be 
interpreted as the ML baroclinic part of the mean velocity. 
The parameterization (Eq. 4a,b) is obtained under condition 
(Eq. lc) and can be obtained from Eq. (7a,b) of CD9 in the 
limit (Eq. lc). 

In Eq. (4a) the velocity Mg may be interpreted as the sub- 
mesoscale induced velocity which is a counterpart of the 
mesoscale induced velocity. As noted earlier, to make the 
analogy with the mesoscale induced velocity more complete 
and since in the mixed layer r, is small due to the strong 
mixing, one may add to Eq. (4a) the term iDg r z , where wf 
is found from the continuity condition: 

3,w+ + V# • m+ = 0 (5) 

The only variable in Eq. (4b) that is not yet parameterized 
is the SM kinetic energy Ksm which we study in the next 


From Eq. (6b), one observes that at the bottom of the ML, 
z=— h, we have that: 

u( — h) — 0, K H (—h) — 0, Fy(—h) — 0 (6c) 

which is a good approximation since SM eddies hardly pen- 
etrate the bottom of the mixed layer (Boccaletti et ah, 2007). 
We also have the additional relations: 

m(0) = 0, k h ( 0) = 0, Fy( 0) = 0 (6d) 

4 MS kinetic energy in terms of resolved fields 

Assuming that the production of Ksm occurs at scales l ~r s 
and since the eddy kinetic energy equation shows that the 
vertical buoyancy flux /-+ acts as the source of Ksm. we em- 
ploy the following relations: 

o 

Ksm — C (r s P K ) 2/3 , P K =<F*> = h~ l J dzF^iz) (7a) 

-h 

In the case of 3-D turbulence where kinetic energy cascades 
from large to small scales and a Kolmogorov spectrum sets 
in, the first relation in Eq. (7a) is simply the statement that 
production=dissipation with the former is defined in the sec- 
ond of Eq. (7a) while dissipation is represented by a Kol- 
mogorov form. In such a case, C=(3/2)Ko(l+A *) where 
A*>0 accounts for the contribution to Ksm of the energy 
spectrum at k<ko and Ko is the Kolmogorov constant whose 
value ranges between 1.4-2. 2. With A*=l, Ko= 2, we have 
C= 6. However, in the 2-D case of interest here, cascade of 
Ksm to smaller scales does not take place. Instead, Ksm 
transforms into SM potential energy which we denote by 
Wsm- Then, in the quasi-stationary case, the production of 
Ksm given in the second of Eq. (7a), approximately equals 
the production of Wsm- Since the latter cascades to smaller 
scales where is ultimately dissipated, we have: 

Wsm^C'^Bk) 273 (7b) 


~ / ~ 
u—y — C7X.11 
Y | 


(4a) 
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where C depends on the spectrum of VV'sm- If we further 
denote by T sm the ratio of the SM kinetic and potential en- 
ergies, from Eq. (7b) we conclude that the first relation of 
Eq. (7a) is satisfied with: 

C = TsmC', Tsm = -Ksm/IEsm (7c) 

It is worth noticing that relations analogous to Eq. (7a,b) hold 
true for mesoscales as well with the proviso that the con- 
stants C , T and therefore C, are different for SM and M 
due to their different dynamics since, as already mentioned, 
SM are trapped in the ML while M form coherent struc- 
tures throughout the entire water column and thus entail the 
dynamics of the deep ocean. Using a mesoscale resolving 
simulation, Eq. (7a) was validated in more than 70 different 
mesoscale resolving simulations (Canuto et al., 2010). In the 
next sections, on the basis of Capet et al. (2008), we estimate 
that C^6. Even though at present we have determined C on 
the basis of only one simulation by C8, we shall show be- 
low that the variable of interest to OGCMs, the tracer flux, 
is only weakly dependent on C. Substituting Eq. (6a,b) and 
rs—Nh/jt |/| into Eq. (7d), we obtain the following alge- 
braic equation for Ksm'- 


5 Wind driven flows 

In this section we study flows driven by strong winds when 
the Ekman velocity exceeds the geostrophic mean velocity 
and compare with results obtained in the submesoscale re- 
solving simulations of Capet et al. (2008). To obtain results 
in an analytical form, we further assume that the ML turbu- 
lent viscosity v~10 -2 m 2 s _1 is z-independent. Under these 
conditions, the mean velocity field can be decomposed into 
geostrophic i/ g and Ekman Up components; with the x axis 
along the wind direction, we have the relations: 

«e = Ae*a(0, v E = (f/\f\)Ae ( l3(0, 

A = (u/)- 1/2 «;, S=z/8e, 8 E = (2v/f) 1 ' 2 

a(f) = cos £ + sinf , f (£) = - 9a(?)/3f (8a) 

where pu 2 is the surface stress and 8 E is the Ekman layer’s 
depth. Below we analyze flows driven by winds of different 
directions with respect to the geostrophic component of the 
mean velocity. 

5.1 Down-front winds 


^SM 2 = C 3/2 (l + P 2 ) _1 hrs(V - Y(f/\f\)e z x V) (7d) 

o 

• V H b, V — h~ 2 J zu(z)dz 

-h 

which is however not convenient for the computation of K$m 
since the latter enters into the right hand side of Eq. (7d) 
through the variables y, as one observes from Eq. (4b). From 
Eqs. (7d) and (4b), one derives the following equation for y: 

A 4 y 4 - A 3 y 3 - y 2 - 1=0, y > 0, N 2 > 0, (7e) 

where: 

A 4 =7r 2 (2C) 3 / 2 (//|/|)( ez xV*). S , (7f) 

A 3 =7r 2 (2C) 3/2 V*-s, F* = V/(h\f\), s=- N~ 2 V H b 


the vector s being the slope of the isopycnals. Equation (7e) 
is valid under the condition 


K _ 2 y 2 K 
Ksm r| f 2 



(7g) 


which is the same as condition (Eq. lc) expressed in terms of 
the resolved fields. We assess this condition in detail in both 
strong and weak wind cases. 

In summary, for OGCMs that resolve M but not SM, the 
parameterization of the z-derivative of the vertical SM tracer 
flux is given by Eqs. (4a, b) and (7e—f). 

To illustrate the solutions we have just derived, in the 
next sections we consider three important cases: (1) strong 
wind driven flows, (2) wind and buoyancy driven flows, 
(3) buoyancy-only driven flows. 


Down-front winds (blowing in the direction of the surface 
geostrophic velocity) drive dense water over buoyant one and 
provide favorable conditions for the generation of subme- 
soscales (Thomas, 2005; Thomas and Lee, 2005). We have: 

«g = mo + SgZ, Vg = 0, mo, Sg > 0 (8b) 

which corresponds to an horizontal buoyancy gradient given 

by: 

V E b = - fS g e y (8c) 

It is worth noticing that S g is a scalar since e y =e z xe x is the 
cross product of the two vectors e z (the unit vector along the 
Earth radius) and e x (the unit vector along the wind direction) 
which is a pseudo-vector. To obtain the submesoscale flux 
and its z-derivative, we need to compute u and u defined in 
Eqs. (4b), (6b), where Ti—up+Ug- Assuming that the mixed 
layer thickness h is considerably larger than the depth of the 
Ekman layer so that the Ekman number E—8 2 E h~ 2 <^\ and 
c'Zo/'A-Cl, from Eq. (8a,b) we derive that: 

0 

<m> = /G 1 / Ti(z)dz, (8d) 

—h 


<U> = U o 


1 

-SJi, <v> 

2 g 


— AE 1/2 

I/I 


u=Ae<a(0+Sg(z+ l -h), u=^A|V/j(C)+£ 1/2 ] (8e) 

A strong wind (larger than the geostrophic wind) is repre- 
sented by the relation: 


A/h » Sg 


(8f) 
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Let us now study the z-derivative of the submesoscale 
buoyancy flux which we obtain substituting Eq. ( 8 c, e) into 
Eq. (4a) with Y—b. The result is: 

3^v= l/IV^l + y 2 )” 1 (9a) 

[e^(0+E 1/2 ]-y e^a(0+A~ l S g (z+^h) j 

As one can check, the z-derivative of this expression is neg- 
ative which implies that the SM vertical buoyancy flux re- 
stratifies the ML. Let us compare the latter effect with that 
of a down-front mean wind that de-stratifies the ML since 
“Ekman flow advects dense water over light” (Thomas et 
ah, 2008). In the approximation of strong mixing, in the 
z-derivative of the mean advection term u-Y7 H b, the largest 
contribution comes from the baroclinic term u-YJ H b (in fact, 
when we differentiate Yi-YJ H b w.r.t. z, the mean velocity 
ii—<u>+u may be substituted by the baroclinic component 
u since the z-derivative of Y/ H b may be neglected). From 
Eq. ( 8 c,e) we have: 

u ■ YJ H b — - |/|S g A[V/l(f) + Zs 1 / 2 ] (9b) 

which we compare with Eq. (9a) in the case of a very strong 
wind: 

Ro 1, i.e. y <$( 1 (9c) 

Under this condition, we may approximate Eq. (9a) as fol- 
lows: 

3 Z F* « \f\S g A[e*M) + £ 1/2 ] = -u-V H b (9d) 

The implication of this relation is that the re-stratification 
by SM largely compensates the de-stratifying effect of the 
mean flow, a conclusion in agreement with Mahadevan et 
al. (2010). To express the eddy kinetic energy in terms of 
the resolved fields, we first find V from Eq. (7d): 

v = y x<?x + v y e y , y x = \ (ae +^s g V ( 10 a) 

Vy = - -A — E l/2 ( 1 - £ 1/2 ) 

y 2 |/| V ) 

Substituting Eq. (10a), together with Eq. ( 8 c), into Eq. (7d), 
we obtain: 

=— (2C) 3/2 (l+r) _1 rs///|S g (10b) 

[(f/\f\)V y -yV x ] 

We recall that in real flows the mixed layer depth exceeds the 
Ekman one and thus in Eqs. (10a) we have £<1. Inspec- 
tion of relations Eqs. (10a,b) shows that both y x y contribute 
with the same sign to the SM kinetic energy. In addition, 
under condition (Eq. 8 b) that S g > 0, these contributions are 

www. ocean-sci .net/6/67 9/20 1 0/ 


positive. Thus, we conclude that the down-front wind gen- 
erates the most vigorous SM eddies, in agreement with the 
results of Thomas (2005), Thomas and Lee (2005), Thomas 
et al. (2008) and Mahadevan et al. (2010) that down-fronts 
winds provide the most favorable condition for SM genera- 
tion. 

As discussed in the previous section, the variable y is so- 
lution of Eq. (7e) in which A 3 4 are given by (making use of 
Eqs. 7f and 10a): 

2A4 = 7t(2C) 3 / 2 4m AE+ -hS g Y 2A 3 (10c) 

N 2 h \ 6 V 

=-^(2C) 3 / 2 ^AE 1 / 2 f 1 -E 1 / 2 ) 

N z h V / 

In the case of a strong down-front wind satisfying condition 
(Eq. 9c) and with £<$(1, we have | A 3 1 > | A 4 1 and the first 
term in Eq. (7c) may be neglected. Then, we obtain: 

Y » (-A 3 )“ 1/3 = (2 C )-'/ 2 E~ 1/6 (10d) 

Next, we discuss condition (Eq. 7g). Using Eq. ( 8 e) for the 
case of a strong down-front wind, we obtain: 

K = -i-A 2 (£ 1/2 - 2£) (lOe) 

Substituting this result into Eq. (7g) together with Eq. (lOd), 
we obtain: 



which is amply satisfied in the simulations of Capet et 
al. (2008). 

5.2 Up-front wind 

Up-front wind (blowing in the direction opposite to that of 
the surface geostrophic velocity). In this case, the contribu- 
tion of the baroclinic component of the mean flow to Eq. (la) 
is given by Eq. (9b) and has the sign opposite to that in the 
previous case since now S g < 0 , i.e., in this case u leads to a 
re-stratification of the ML. 

It is worth treating strong winds first. In this case, the first 
term of V x in Eq. (10a) dominates, we have y x >0, Vy<0 and 
thus the expression in the square bracket of Eq. (10b) is neg- 
ative. Using the relation (l+y 2 ) -1 = Ksm/{Ksm + rg f 2 ), 
we conclude that the only solution of Eq. (10b) is Esm= 0 
which implies that submesoscales are not generated. 

On the other hand, in the case of weak winds when sub- 
mesoscales are generated, they re-stratify the ML. Thus, the 
re-stratification effect of the mean flow is strengthened by 
that due to SM, a conclusion in agreement with the results of 
Mahadvan et al. (2010). 
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5.3 Wind perpendicular to the geostrophic velocity 

In such flows we have: 

Vg — vo + SgZ, Ug — 0 (11a) 

and thus: 

V fl £=|/IVx (Hb) 

In contrast with Eqs. (8d-e) and (10a), in this case, the 
geostrophic component contributes to the y-projections of 
<«>, h, and V. In particular, we have 

2V X =AE, 2V y =-|yy AE 1/2 (l-E l/2 ^) + ^hSg (11c) 

Substituting this relation into Eq. (7d), we obtain: 

*S 3 m=C 3/2 (l+K 2 ) ^ _1 rsh \f\S g [V x +y (// 1/|) V y ] (1 Id) 

Inspection of relations Eq. (llc.d) shows that the compo- 
nents Vx y contribute to Ksm with the opposite sign. Still, 
in the case of a strong wind corresponding to a small y, the 
first term in Eq. (lid) dominates which yields a positive K^m 
if S g >0. In this case, the direction of the wind coincides with 
that of the horizontal buoyancy gradient (up-gradient wind). 
Under the same condition Sg> 0, the baroclinic component of 
u de-stratifies the ML since in this case instead of Eq. (9b), 
we have: 

u-VHb=\f\S g Ae ( aU) (lie) 

which has a positive z-derivative and thus de-stratifies the 
mixed layer, while SM tend to re-stratify the mixed layer. As 
one can see from Eqs. (9b) and (lie), in both cases (Sect. 5.1) 
and (Sect. 5.3) “Ekman flow advects dense water over light” 
(Thomas et al., 2008). On the contrary, if V# /? has a direction 
opposite to that of the wind (down-gradient wind), so that 
Sg<0, strong winds tend to re-stratify the mixed layer and do 
not generate submesoscale eddies. 

In conclusion, our analytical results for flows driven by 
wind and baroclinic instability with different directions of 
the wind and the buoyancy gradient, are in agreement with 
the results of eddy resolving simulations by Mahadevan et 
al. (2010). 

6 Testing the SM parameterization against simulations 
with wind and buoyancy 

In this section, we compare our model results with Capet et 
al. (2008, C8). We begin with Fig. 12 of C8 which gives 
d- F 2 ( z ) which we compare with our model (Eq. 4a). To do 
so, we assume that the direction of V# T coincides with that 
of Vtfb and that the wind blows in the down-front direction. 
Then, in analogy with Eq. (9a), from Eqs. (4a) and (8a) we 
obtain: 

a.F^Ad+rr 1 ( 12 a) 


[e f i8(f)+£ 1/2 ]-y^ f a(f)+A- 1 Sg( z +i/i jv ff |r| 

Since from Fig. 11 of C8 we have 
0.8-HT 5 °Cm -1 <| V# T| <2.7-10 -5 °Cm“ 1 , we use 

| V# r| = 1.5-10 -5 °Cm _1 . As for the buoyancy frequency 
N, we use the characteristic value /V=10 -3 s -1 . Next, we 
need the values of A, 5 g , /;, &e and /fsM defined in the 
previous section. From the C8 data presented in Fig. 10, we 
derive the following results: 

S g «6 x 10“ 4 s“', A«6 x 10 -2 ms -1 , (12b) 

&E ~ 10 m, /?~40 m, AIsm ^ 1.8 x 10 -3 m 2 s -2 

While for the determination of A and Sg we need the profiles 
of the projections of the mean velocity u(z ), in their Fig. 10, 
C8 present only the absolute value \u(z)\. However, assum- 
ing that the wind blows in the down-front direction, one can 
extract A and S g from their Fig. 10. Finally, from Eqs. (12b), 
(4b) with rs=(N/n\f\)h and Eq. (lOe), we derive that: 

y = 0.3, K = 1.1 x 10“ 4 m 2 s“ 2 , (12c) 

E = ( S E /h ) 2 = 1/16 

Using this value of y in Eq. (7e) together with Eq. (10c), we 
derive that: 

C = 6 (12d) 

Regrettably, the C8 data are the only ones available that al- 
lowed us to estimate the parameter C and other data would 
be welcome to confirm Eq. (12d). Fortunately, as relation 
Eq. (lOd) shows, the dependence of y on C is rather weak. 
In addition, under condition (Eq. 9c), which is satisfied by 
the first result (Eq. 12c), in Eq. (12a,e) the terms linear in y 
are small compared to the main terms which do not contain 
y . Thus, we conclude that in the case of a strong down-front 
wind, the effect of C on the SM tracer flux is relatively weak. 
Next, substituting Eq. (12b,c) into Eq. (lc), one can see that 
the condition is amply satisfied. Finally, in Fig. 1 we com- 
pare the z-profile of —d-F 2 (z) from Eq. (12a) with that of 
Fig. 12 of C8 (dashed line). In Fig. 2, we compare the pro- 
files of the fluxes F 2 (z) from the present model: 

^v=(i+y 2 ) _1 A ( 12e ) 

cost;+2z/ h)&E~y &Ee i; sin/+^A _1 S' g z(z-zo) j 

\v h t\ 

with that of C8 which we compute using C8 data for the z- 
derivative of the flux shown in Fig. 12. As one can observe, 
the profiles of the fluxes are quite close throughout the mixed 
layer depth. As for the profile of —d z Fy (z), they are quite 
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Fig. 1. —d z Fy(z) for ICCO simulation data by Capet et al. (2008) 
(dashed curve taken from their Fig. 12) and for the present model 
(Eq. 12a) (solid curve). 



10"Fj(°Cms“‘) 

Fig. 2. Profile of the vertical heat flux Fy(z) for the ICCO simu- 
lation data by Capet et al. (2008) (dots) obtained by integrating the 
dashed curve in Fig. 1, and for the present model (Eq. 12c) (solid 
curve). 



Fig. 3. Absolute value of the mean velocity considered in the 
present analysis which corresponds to the case of a down-front wind 
and which is given by Eqs. (8a, b), (12b) (solid curve) and the one 
computed in ICCO simulation of Capet et al. (2008) (dashed line 
which is taken from their Fig. 10). 

no wind (Fox-Kemper et al., 2008, FFH) which can also be 
used to test of our model predictions. Following FFH, we as- 
sume that the mean velocity is in thermal wind balance with 
the mean buoyancy field and that the mean buoyancy gradient 
does not vary inside the mixed layer, i.e., u z =f~ l e z x 
is z independent. Irrespectively of the surface value of ii, 
from Eq. (4b) and the second of Eq. (6b), we derive that: 

~ 1 _ 1 

u= — (2z+h)e z xV H b, u=—(z+h)e z xV H b (13a) 

Taking Y—b in Eq. (6a,b), and using Eq. (13a), the buoyancy 
flux is given by: 

F v (z) = <l>(y)-^— (l -£ 2 )|V^| 2 , (13b) 

16|/| V / 


close in the upper half of the mixed layer but differ some- 
what in the lower half. We think this is due to the similarity 
of the mean velocity profile (Eq. 8a) with that in the C8 at 
small depths and by an unavoidable difference in the lower 
part of the mixed layer due to the different profiles of the ver- 
tical viscosity used here and in the C8 simulations. While in 
our analysis we adopted an Ekman profile corresponding to 
a v(z) = const., C8 adopted a more realistic model for v(z). 
In spite of that difference, the profiles |m(z)| compare well, 
as seen in Fig. 3. 


1 = 1 + 2 zh-\ <t>(y)= 

1 + y 

with: 

F v (0)=0, F v (-/0=0, Fy(-h/2)=<i>\V H b\ 2 >0 (13c) 

To complete this parameterization of Fy, y must be obtained 
from solving Eq. (7e). Substituting Eq. (13a) into the second 
of Eq. (7d), we obtain: 

h 

V = — — e z x V/fb (14a) 


7 No-wind case - Buoyancy only 


Substituting this result into Eq. (7f), we obtain: 


In addition to the realistic case studied by C8, there are data 
from simulations corresponding to the less realistic case of 


A 4 = — (2C) 3/2 Ri~ 1 , A 3 = 0, Ri = ^ N 
12 \V H b\ 


(14b) 
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Fig. 4. Richardson number dependence of the two ratios in 
Eq. (15a, d). The first ratio, represented by full circles with error 
bars, correspond to updated results from Fox-KempeT Comments 
to our manuscript during the Discussion phase: Ocean. Sci. Dis- 
cuss. 6, C916-C926, 2009, http://www.ocean-sci-discuss.net/6/. 
The results of the present model given by Eq.(15e) with C=6 are 
shown by the full line. The blue line corresponds to the FFH pa- 
rameterization, <t>=l. 


The solution of Eq. (7f) is then as follows: 

y~ = C* ^Ri + jRi 2 + 2 Ri/C^j , C* = ^ (2C) 3/2 (14c) 

If we employ the value C=12 that we have determined from 
the data of Capet et al. (2008), we can compute the function 
(Eq. 13b): 

<S>(y)=Q(Ri) (14d) 


which we exhibit in Fig. 4. We further notice that in spite 
of the fact noted before that we have only one set of data to 
determine C, Fig. 5 shows that different C’s have only an 
overall marginal effect on the buoyancy flux in the interval 
l<Ri<100. 

Next, we find the limits of applicability of the parameter- 
izations (Eqs. 13b and 14c, d). To this end, we find the baro- 
clinic mean kinetic energy averaged over the ML depth K 
using the definition (Eqs. 7g, 4b and 13a): 

K = ^h 2 f- 2 \V H b\ 2 = (14e) 


This result shows that at the surface the mean baroclinic ki- 
netic energy K is twelve times smaller than the mean kinetic 
energy Kyi—hu 2 . To compute the SM kinetic energy, we use 
the first definition (Eqs. 4b and 14c) and derive that: 

2 

Y = Y(Ri) (14f) 


(hN\ 
2^sm = ( — ) 


\7ty / 



Fig. 5. Same as Fig. 4 but with two different values of the parameter 
C=3 (dashed), 6 (solid) to show the weak dependence on such a 
parameter. 


Thus, condition (Eq. lc) becomes: 
r 2 


Ri > 


JT- 


36(3C 3 / 2 -!) 


- 2 x 10“ 3 


(14g) 


Next, we compare Eq. (13b) with the FFH data and recall 
that in their Fig. 14e the authors plot the ratio: 


A (data) = 


TV (data) 
E v (FFH) 


(15a) 


where: 


Fv(FFH) = 0.06|/r 1 /t 2 /x(z)|V///7| 2 , (15b) 

/*(z) = ( l~f 2 )(l+5§ 2 /21) (15c) 


is the parameterization suggested by FFH by fitting their sim- 
ulation data. To compare the model results with the simula- 
tion data, we construct the ratio: 


A (model) = 


Tv (model) 
TV (FFH) 


(15d) 


Thus, A(FFH)=1, by definition (blue line in Fig. 4). To com- 
pute Eq. (15d) in our model, we notice that the profile p(z) 
in Eq. (15c) has the additional factor (l+5£ 2 /21) in com- 
parison with the profile (Eq. 13b). The difference does not 
exceed 25% and we neglect it when substituting Eqs. (13b) 
and (15b, c) into Eq. (15d). As a result, we obtain: 


A (present model) = <t>( Ri ) 


(15e) 


which we show in Fig. 4 (red line). The present parame- 
terization yields a good representation of the FFH simula- 
tion data especially the curvature exhibited by the data. The 
agreement is somewhat worse at Ri> 10 3 which is due to the 
use of Killworth’s (2005) approximation to neglect r in 
fact, such an approximation becomes questionable when Ri 
is large, that is, when b z =N 2 is large. Finally, we discuss 
whether the FFH flux (Eq. 15b,c) without wind can represent 
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the case with wind. To that end, we take the value of Sg in 
the first of Eq. (12b) as determined from C8 simulations and 
substitute it in Eq. (lib) with the result: 

V H b « 0.5 ■ 10“ 7 s“ 2 (15f) 

Substituting this result in Eq. (15b,c), and using the same 
mixed layer depth h= 40 m as in C8, we obtain: 

max Ev(FFH) = 2.4 • 10“ 9 m 2 s“ 3 (15g) 

If one compares this value with the value from C8 (Fig. 2): 

max Fy — ga^max Fy = 2 • 10 -8 m 2 s -3 (15h) 

one concludes that the FFH flux formula with no wind under- 
estimates the flux by about an order of magnitude, a conclu- 
sion in agreement with Mahadevan and Tandon (2006) who 
stressed that “winds play a crucial role in inducing subme- 
soscale structure” not to mention the multifaceted and im- 
portant implications on the mixed layer stratification caused 
by the down-wind, up-wind vs. buoyancy topology described 
in Sect. 6. 

8 Conclusions 

Recently, there has been a considerable interest in sub- 
mesoscales which are oceanic structures of O (1 km) size 
and life times of the order of days. High resolution numer- 
ical simulations have been the best source of information 
to assess the parameterizations of SM fluxes to be used in 
OGCMs that do not resolve such features. If one consid- 
ers that the highest resolution of about 1/10° in stand-alone 
OGCMs can represent structures of about 10km size which is 
10 times larger than the SM sizes and that OGCMs employed 
in thousand years runs for climate studies have resolution of 
about 1° (corresponding to structures 100 times larger than 
SM), it seems clear that a good deal of important physical 
processes have thus far gone unrepresented in those OGCMs. 

The present work presents an analytical parameterization 
of the SM vertical flux of an arbitrary tracer. The main fea- 
tures can be described as follows: 

(a) no other SM parameterization exists (to the best of our 
knowledge) that provides an analytical expression for 
the vertical flux of an arbitrary tracer under arbitrary 
buoyancy and wind (strength and direction), 

(b) the SM parameterization by FFH does not include winds 
and it is limited to the buoyancy field which cannot be 
used to describe tracers such as the ones needed in the 
C-cycle models in OGCMs for climate studies, 

(c) the results of the realistic simulations by C8 and MTF 
with baroclinic instabilities and winds, are well repro- 
duced by our model. 


(d) our parameterization given by Eqs. (4a-c), (7e-f), 
(12d), is to be used in OGCMs that resolve mesoscales 
but not sub-mesoscales, 

(e) in a different study (Canuto and Dubovikov, 2010), we 
have derived the parameterization for the tracer verti- 
cal flux for OGCMs that do not resolve either subme- 
soscales or mesoscales. 

Appendix A 

The non-linear terms g H 

As discussed in textbooks on Turbulence theories (e.g., 
Batchelor, 1970; Fesieur, 1990; McComb, 1992), the 
stochastic Fangevin equation has played a major role in tur- 
bulence modeling studies (Kraichnan, 1971; Feith, 1971; 
Herring and Kraichnan, 1971; Chasnov, 1991). Though most 
turbulence models are presented in terms of the energy spec- 
trum, which is a second-order moment, the starting point is 
always the Navier-Stokes equations (NSE) presented in the 
form of a stochastic Fangevin equation in k-space: 

d t u\(kj) = - v d (k)k 2 Ul {k,t) + f^(kj) (Al) 

in which the non-linear (NL) term of NSE is represented by 
the two terms: the turbulent forcing f*(k,t) which is due to 
the infrared (small k ) part of the NF interactions and ultravi- 
olet part which is represented by the enhanced ^-dependent 
dynamical viscosity v d (k)—v+v t (k), where v is the kine- 
matic viscosity while v t (k) is a turbulent viscosity discussed 
in Appendix B. As discussed in the references cited above, 
the dynamic equation for the energy spectrum E(k) is ob- 
tained by multiplying Eq. (Al) by u* ik') and integrating over 
n=k/|k|, the result being: 

d,E(k) = Mk ) - 2 k 2 v d (k)E(k) + A ex , (A2) 

where the work A t is given by: 

(A3, 

On the other hand, the general equation for E (k) is given by 
(Batchelor, 1970, Eq. 6.6.1) 

d,E(k) = T(k) - 2 vk 2 E(k) + A ext (A4) 

where T (k) is the non-linear transfer. From Eqs. (A3.A4) it 
follows that: 

T(k) = A { (k) -2v t (k)k 2 E(k) (A5) 

Within the closure model developed by Canuto and 
Dubovikov (1997), the form of A t (k) is given by: 

k 

9 E(k) i r , 

A t (k) = -r(k) , -r(k) — J p~v,(p)dp (A6) 

o 
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The key feature of this closure is that A { (k) is proportional 
to the derivative d^E which vanishes in the vicinity of the 
wavenumber k—ko where E (k) has its maximum. This re- 
duces the two NL terms in Eq. (Al) to the second one only 
which, in the notation of Eq. (2d,e), implies that: 

Qh = *%v d (ko)u', v d = xd (A7) 

Appendix B 

Turbulent viscosity v t (k) 

Contrary to the kinematic viscosity v which does not depend 
on the size of the eddies, the turbulent viscosity v, ( k ) which 
is due to the NL interactions, depends on the eddy size and its 
sum to v is called the dynamical viscosity, v d (k)—v+v t (k). 
The search for a suitable expression for v t (k) dates back 
many decades and the first explicit expression is the heuris- 
tic one proposed by Heisenberg as discussed in Batchelor’s 
book (1970, Sect. 6.6, Eq. 6.6.13), 


which is the well-known Richardson 4/3 law diffusivity 
~£ 4 / 3 . Finally, relation (Eq. B2) shows that there is no such 
a thing as a unique turbulent viscosity since each eddy feels 
its own turbulent viscosity. In Eqs. (2e) and (3a) we are inter- 
ested in the function v t (k)^v d (k) in the vicinity of the max- 
imum of the energy spectrum k—k^. Assuming that most of 
the energy is contained in that region, from Eq. (B2) we de- 
rive that v d ~k^ 1 K^ 2 . Thus, from Eq.!(A7) it follows that: 

Qh = *£«' (B5) 

which is the closure form in Eq. (3a). The closure for the 
tracer field is analogous. 

Appendix C 
Derivation of Eq. (4) 

Cl Sub-mesoscale tracer field x' 


OO 

v,{k) = Y J p- y2 E(p)^ 2 dp y= 0(1) (Bl) 

k 

where E (k) is given by Eq. (A2) is the kinetic energy spec- 
trum whose integral over all wavenumbers yields the eddy 
kinetic energy Ke . As discussed by Batchelor, Eq. (Bl) was 
successfully used to derive the Kolmogorov spectrum. A non 
heuristic derivation of v, (k) has however been lacking un- 
til recently with the advent of methods to treat the Navier- 
Stokes equation borrowed to a large extent from quantum 
field theory. A full presentation was made by the present 
authors in 1997 with the final result: 

/ oo \ 1/2 

V t {k) = y + \j P ~ 2 E(p)dp 1 (B2) 

Equation (B2) has several interesting features worth dis- 
cussing. First, it says that an eddy of arbitrary size (~k _1 ) 
feels the effects of all the eddies smaller than itself, as the in- 
tegral begins at k and accounts for all the wavenumbers from 
k to infinity. Equation (B2) naturally reduces to the kine- 
matic viscosity v when the size of the eddy is very small and 
k tends to infinity. Due to the presence of the kinematic vis- 
cosity, Eq. (B2) is valid for arbitrary Reynolds number since 
it can be rewritten as: 

r 9 1 1/2 fK 1 / 2 

v,(k) = v\l + Re(k) 2 J , Re{t) — (B3) 

If one employs the Kolmogorov spectrum 
E(p)—Kos 2 ^k~ 5 ^, one obtains in the large Re regime: 

( t 3Ko e 2 / 3 \ 1/2 ,,, in 

Re » 1 : v,(k) = I v 2 + — — j « g 1 / 3 ^ 4 / 3 (B4) 


Substituting Eq. (3a) into the tracer Eq. (2e), we obtain the 
following expression for the submesoscale tracer field: 


«' ■ Vflt 
X +i(k ■ Ti — co) 


(Cl) 


where |k|=ko=r^“ 1 . As in the case of mesoscales discussed 
in CD5, the frequency co is obtained by solving the eigen- 
value problem mentioned below Eq. (3a) with the following 
result (dispersion relation): 


co(k) = k ■ u d (C2) 

This relation can be interpreted as the Doppler transforma- 
tion for the frequency provided that in the system of coor- 
dinates moving with velocity u d , the submesoscale flow is 
stationary, in which case co— 0. Stated differently, relation 
(Eq. C2) implies that u,j is the eddy drift velocity whose ex- 
pression in terms of mean fields is analogous to that for the 
case of mesoscales given in Eq. (4f) of CD6: 

1 9 

u d = <u> + -€“e z x (/? — / <d z s>) (C3) 

where s—— Vffb/N 2 is the slope of isopycnal surfaces and 
ft=V f . The bracket averaging is defined as follows: 

o o 

<•> = f (.)K^ 2 (z)Jz/ J K ] H /2 (z)dz (C4) 

—h —h 

Due to the smallness of the scale i characterizing subme- 
soscales, the second term in Eq. (C3) is negligible and thus 
in good approximation we have that: 


u d — <u> 


(C5) 
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which changes Eq. (Cl) to the form: 


r 


/ 


u' ■ V# r _ _ 

— — , u = u — <«>, 

X +ik ■ u 


x = r'K? 


(C6) 


Relation (Eq. C2) implies that the dependence of the subme- 
soscale fields on co is of the form: 


A'(co,k) = A'(k)8(w-k-u d ) (C7) 

and therefore in the (f, &)-space the fields A ' depend on time 
as follows: 

A f (t,k) = A'(k)exp(— ik ■ ujt) (C8) 

Due to relations (Eqs. C7.C8), after substituting Eq. (3a) in 
Eq. (2e), the latter may be solved in both (&>, k) and (t, k) 
representations. 


C2 Sub-mesoscale velocity field w' 


It is convenient to begin by splitting the mesoscale velocity 
field ii into a rotational (divergence free, solenoidal) and a 
divergent (curl free, potential) components: 

u'(k) — UR(k) + ui)(k): u R (k) — n x e z u R (k), (C9) 


iiD(k) — nt/D, n — k/k 


and thus the third equation in Eq. (2e) becomes: 

3 z w’ — — ik ■ u = —ikiiD (CIO) 


To determine mr.d. we substitute the second relation (Eq. 3a) 
in the second equation in Eq. (2e) and derive the following 
expressions: 

ud = f~ 1 (x + ik -u)u R (Cl 1) 


MR = 


ikp 1 p' 


1 + / ~(X + ik ■ u) 


U — u 


<u> 


(Cl 2) 


These relations, as well as Eq. (C6), are valid in both (&>, k) 
and (r , k) representations. Below, we will use them in the 
(f,k) representation together with Eq. (C8). Due to the third 
relation of Eq. (3a), the further use of Eqs. (Cll) and (C12) 
is considerably simplified under the condition 

Ksm/K » 1 (Cl 3) 


which allows us to neglect the second terms in the brackets 
in Eqs. (Cll) and (C12). Condition (Eq. Cl 3) coincides with 
Eq. (lc). 


C3 Spectrum of the vertical SM flux 

In the dynamical Eq. (la) one needs d z Fy which we shall 
write as follows: 

3 z Fy — w' z t' + w' x ’ z ~ w' z r' (C14) 

where in the last expression we have neglected the term w' r' z 
in accordance with the adopted approximation r' and be- 
cause it is of a higher order in z. Substituting Eq. (Cll) into 
Eq. (CIO), we obtain the expression for 3 z w' which allows 
us to compute 3 z Fy using Eq. (C14). 

The strategy of computing submesoscale fluxes, which are 
bilinear correlation functions, consists in computing these 
functions in (t,k)- space which, in the approximation of ho- 
mogeneous and stationary mean flow, have the form: 

A'(t,k')B'*(t,k) = A' B'*(k)S{k - k') (C15) 

which, because of relation (Eq. C8), does not depend on t. 
The function Re A’ B'*(k) is usually referred to as the density 
of A'B' in k-space. The spectrum of the correlation function 
A!B' is: 

~MW{k) — [ ReA'B'*(k)8(k - \k\)d 2 k (C16) 


i.e., the spectrum of A'B' is obtained by averaging 
Re A' B'*(k) over the directions of k and multiplying the re- 
sult by jtk. Finally, the correlation function A'B' in physi- 
cal space is obtained by integrating over the spectrum. Fol- 
lowing this procedure, from the first of Eq. (C6) and using 
Eqs. (C9.C10), we derive the relation: 

Re w' z x'*(k) = (C17) 


- Im Jk/ 2 (x + fA:-M)|^MDMR(A:)nxe z -|-|MDl 2 (*)n]}- V#r 
Next, from Eqs. (Cl 1,02) we derive the following relations: 
Re^l(k) = X r l JuRf(.k), (Cl 8) 

Im uou^ik) = f~ l k ■ m|mr| 2 (A:) 

JuD^(k) = x 2 r 2 V^-(k), Wik) (C19) 

= |wdI 2 (AO + |mr| 2 (AO = (1 + x 2 / _2 )|mrI 2 (£ ) 

Substituting Eq. (Cl 8) into Eq. (C17) and averaging over the 
directions of k, we obtain the spectrum of w' z x '(k). Under 
conditions (Eq. Cl 3), we get: 

w'lX’ik) (C20) 

= - £oX _ 2 [x/ _ 1 |mr| 2 (A:)h xe + \u D \ 2 (k)u^ - V H x 
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where: 

nk\u R \ 2 (k) = 2(l + X 2 f~ 2 ) * E(k), 7tk\u D \ 2 (k) (C21) 

= 2(l + X- 2 / 2 )“%), E(k) = \itk\^(k) 

where E(k) is the spectrum of the total (rotational + diver- 
gent) eddy kinetic energy. Due to relation (Eq. C14), the 
left hand side of Eq. (C20) multiplied by jvk is the spectrum 
of the z-derivative of the vertical flux 3 z F\(k). Thus, mul- 
tiplying Eq. (C20) by nk, using Eq. (C21), we obtain the 
following expression for the spectrum of 3 z Fy(k) near the 
maximum of the energy spectrum: 

3 z Fy(k) = *(*)• Vtf r (C22) 

where: 

4>(k) = - 2E{k)klx~ 2 (C23) 

[x/ _1 d + X 2 /” 2 )” 1 " x e z + (1 T X _2 / 2 r'«] 

Assuming that the spectra Fy(k) and E(k) have similar 
shapes, integration of Eqs. (C22,C23) over k reduces to the 
substitution of E(k) and F\(k) with the eddy kinetic energy 
A"sm and F\ in physical space. Thus, we get Eq. (4a,b). 
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